The Kelvin–Voigt equation in isotropic viscoelastic solids has been successfully employed to simulate the tcFUS field with mode conversions in the skull.^[32,34] In an elastic material, the strain at each point is only dependent on the instantaneous local stress. The stress and strain are related by the stiffness. For an anisotropic medium, this relationship can be written using Einstein summation notation as

To describe the propagation of elastic waves, the momentum conservation relationship between stress and strain must be combined, which is written by the stress and the particle velocity as

Equation (4) can similarly be written as a function of the stress and particle velocity using Eq. (2) as follows:

Equations (7) and (8) constitute the Kelvin–Voigt equation and are used to calculate the tcFUS field.

The wave field separation is based on the fact that a compressional wave has no vortex and a shear wave has no divergence either. The particle velocity is separated into contributions of compressional and shear waves denoted respectively as and .

Since the compressional and shear parameters of the skull are different, the heat deposition in the skull is contributed by the separated compressional and shear waves above. For single-frequency ultrasound used in tcFUS, the heat deposition can be written as

The Pennes bio-heat equation^[37] is commonly used to describe the temperature field and can be written as

The staggered-grid finite difference time domain (FDTD) method^[38,39] is employed to numerically solve the Kelvin–Voigt equation and the wave field separation equations. The parameter averaging technique^[40] is used to satisfy the continuity conditions on the interfaces of solids (skull) and fluids (water and brain inside the skull). The non-splitting convolutional perfectly matched layer (NCPML)^[41,42] is utilized to eliminate the numerical reflections at the numerical boundary. The FDTD method is also used to numerically solve the Pennes equation.

Generally, computed tomography (CT) scans of a monkey skull are used to reconstruct the three-dimensional (3D) skull model, as shown in Fig. 1. There are totally 264 CT scans and each CT scan has 512×512 pixels. The original CT scans have 512×512×264 voxels with a voxel size of 0.3 mm×0.3 mm×0.3 mm, and the linear interpolation is needed in most of the situations to satisfy the numerical requirement of the FDTD. To obtain the strongly inhomogeneous acoustic parameters of the skull, the relationship between the acoustic parameters and CT values of the skull is adopted as follows:^[43]

Fig. 1.

	Figure Option ViewDownloadNew Window
	Fig. 1. Reconstruction of three-dimensional skull model.

The planar phased array for tcFUS is shown in Fig. 2(a). The array is comprised of square elements, each with a size of 3.2 mm×3.2 mm and a kerf width of 0.3 mm. The designed focus is located at 60 mm in front of the array center. The frequency of the array is 0.8 MHz. For each element of the array, the acoustic pressure radiated is

Fig. 2.

	Figure Option ViewDownloadNew Window
	Fig. 2. (color online) Schematic diagrams of (a) planar phased array, and (b) cross section across central axis with focus.

Eight computational models are chosen as shown in Fig. 3. In Figs. 3(a)–3(e), θ is 0° and FD increases from 10 mm to 50 mm in steps of 10 mm, which indicates the relative location changes of the skull on the central axis. Fixing FD at 30 mm, we rotate the skull clockwise around the rotation point in Fig. 2(b) by 10°, 20°, and 30° as illustrated in Figs. 3(f)–3(h), respectively.

Fig. 3.

	Figure Option ViewDownloadNew Window
	Fig. 3. (color online) Computational models where ((a)–(e)) θ = 0° and FD increases from 10 mm to 50 mm in steps of 10 mm, and ((f)–(h)) FD = 30 mm and θ increases from 10° to 30° in steps of 10°.

The time-reversal method^[44] is used to precisely focus ultrasound from the array through the skull at the focus. Ultrasound emitted from a virtual sinusoidal point source at the focus is recorded by each element of the array. The signal received by an arbitrary element is chosen as the reference set as and then undergoes self-correlation and cross-correlations with the received signals of any other element to calculate the initial emitting phases as follows:

The initial temperature is . The acoustic and thermal parameters used in the simulation are listed in Table 1. By taking Fig. 3(c) for example, the distributions of the CT values and the acoustic and thermal parameters are shown in Figs. 4(a)–4(h).

Fig. 4.

	Figure Option ViewDownloadNew Window
	Fig. 4. (color online) Snapshots on Fig. 3(c) for (a) CT values, (b) density, (c) compressional wave velocity, (d) shear wave velocity, (e) compressional absorption coefficient, (f) shear absorption coefficient, (g) specific heat, and (h) thermal conductivity.

Table 1.

Table 1.

Table 1. Acoustic and thermal parameters used in the simulation. . Sound Absorption coefficient Density Thermal conductivity Specific velocity/(m/s) /(Np/m, 1 MHz) /(kg/m3) /(W/(m/K)) heat/(J/(kg/K)) 1500 7.37 1000 0.54 4180 3100 293.12 2200 1.30 1840 1500 0.000288 1000 0.52 3700 1560 3.91 1030

Table 1. Acoustic and thermal parameters used in the simulation. .

What needs to be pointed out is that complex guided waves may be excited in the skull, because the thickness of the skull is close to the wave length, multiple reflections and refractions occur at the boundaries of the skull and the inhomogeneity of the skull leads to strong mode conversions and coupling between compressional and shear waves in the skull. Nevertheless, the guided wave is another form of compressional wave or shear wave in the skull, which has been included in the Kelvin–Voigt equation as the Kelvin–Voigt equation can describe the full wave field in isotropic solids. The guided waves in the skull may make sense in monitoring the tcFUS therapy and will not be discussed here since what we are going to study in this paper is not related with the guided waves.

The normalized acoustic intensity ( , ) distributions without and with mode conversions on the cross sections in Fig. 3 are shown in Figs. 5 and 6, respectively. The and denote the acoustic intensities contributed by compressional and shear wave, respectively. The is the acoustic intensity radiated by the array. To obtain the results without mode conversions in Fig. 5, and are both set to be 0. Besides, the normalized acoustic intensity distributions without and with mode conversions on the central axis are drawn from Figs. 5 and 6 for contrast as plotted in Fig. 7.

Fig. 5.

	Figure Option ViewDownloadNew Window
	Fig. 5. (color online) Normalized acoustic intensity distributions on the cross section without mode conversions where ((a)–(e)) θ = 0° and FD increases from 10 mm to 50 mm in steps of 10 mm, and ((f)–(h)) FD = 30 mm and θ increases from 10° to 30° in steps of 10°.

Fig. 6.

	Figure Option ViewDownloadNew Window
	Fig. 6. (color online) Normalized acoustic intensity distributions on the cross section with mode conversions, where ((a)–(e)) θ = 0° and FD increases from 10 mm to 50 mm in steps of 10 mm, and ((f)–(h)) FD=30 mm and θ increases from 10° to 30° insteps of 10°.

Fig. 7.

	Figure Option ViewDownloadNew Window
	Fig. 7. (color online) Normalized acoustic intensities on the central axis with (red lines) and without (green lines) mode conversions, where ((a)–(e)) θ = 0° and FD increases from 10 mm to 50 mm in steps of 10 mm, and ((f)–(h)) FD=30 mm and θ increases from 10° to 30° in steps of 10°.

From Figs. 5–7, the focal regions for all computational models generally lie around the designed focus. By extracting the coordinates where the maximal normalized acoustic intensity in the focal region (namely the focusing gain) is located, the focusing deviations between the real and the designed focuses are given in Fig. 8(a). The average focusing deviations with and without mode conversions are 0.3187 mm and 0.3629 mm, respectively, indicating that the time-reversal method can guarantee the precise focusing of the transcranial ultrasound. Nevertheless, the focusing deviation tends to increase with the increase of the rotation angle of the skull. As can be seen from Figs. 5(f)–5(h) and Figs. 6(f)–6(h), when θ increases, the shape of the focal region is stretched towards the rotation direction of the skull, especially when θ = 30°, which may account for the increase of the focusing deviation.

Fig. 8.

Figure Option ViewDownloadNew Window

Fig. 8. (color online) Y: with mode conversions; N: without mode conversions; NAI: normalized acoustic intensity. Computational model numbers 1–8 represent the computational models in Figs. 3(a)–3(h), respectively. (a) Focusing deviation between real and designed focuses, (b) focusing gain of the acoustic intensity in the focal region and the maximal normalized acoustic intensity around the skull, (c) drop percentage of the focusing gain with mode conversions compared with focusing gain without mode conversions, and (d) ratio of the maximal normalized acoustic intensity around skull to focusing gain.

According to Figs. 5–7, the focusing gains of the acoustic intensities in the focal regions are given in Fig. 8(b). In either computational model, the focusing gain in the focal region with mode conversions is lower than that without mode conversions. The drop percentages of the focusing gains with mode conversions compared with the focusing gains without mode conversions are given in Fig. 8(c) and the average drop percentage for all computational models is 12.05%. For θ = 0° and FD = 10 mm–50 mm, the focusing gain in Fig. 8(b) decreases with a slow tendency. The focusing gains with and without mode conversions are generally around 30.0 and 33.5, respectively. The drop percentage in Fig. 8(c) firstly decreases slightly and subsequently increases slowly, with a 10.8% average value. The fluctuation of the differences between the focusing gains with and without mode conversions is small. For FD = 30 mm and θ = 10°–30°, the focusing gain declines rapidly, from 30.9 to 8.1 without mode conversions and from 27.5 to 6.9 with mode conversions. The difference between the focusing gains with and without mode conversions becomes smaller and smaller. The drop percentage reaches a maximum of 16.7% when θ = 20° and an average value of 14.15% for FD = 30 mm and θ = 10°–30° is greater than 10.8% for θ = 0° and FD = 10 mm–50 mm. When θ increases, neglecting mode conversions will cause more acoustic energy to be reflected in comparison with considering mode conversions, which will result in the smaller difference between the focusing gains with and without mode conversions. We believe that when θ continues to increase, the focusing gain with mode conversions will exceed that without mode conversions, leading to a negative drop percentage.

Figures 5(a)–5(e) and figures 6(a)–6(e) show that the stronger normalized acoustic intensities are mainly concentrated around the middle and both sides of the skull for θ = 0° and FD = 10 mm–50 mm, of which the maximal ones are located around the middle of the skull. For FD = 30 mm and θ = 10°–30°, although the normalized acoustic intensities around the middle of the skull are also strong as shown in Figs. 5(f)–5(h), Figs. 6(f)–6(h) and Figs. 7(f)–7(h), the maximal ones are focused near the upper side of the skull which is rotated towards the focus. The maximal normalized acoustic intensities around the skull for all computational models are also picked up and plotted in Fig. 8(b). For θ = 0° and FD = 10 mm–50 mm, the maximal normalized acoustic intensities around the skull with mode conversions are higher than those without mode conversions and the differences between them fluctuate less. Besides, the maximal normalized acoustic intensities around the skull are lower than the focusing gains. When FD = 10 mm, the maximal normalized acoustic intensities around the skull with and without mode conversions, especially with mode conversions, are close to the focusing gains. A sharp drop between the maximal normalized acoustic intensity around the skull for FD = 10 mm and the one for FD = 20 mm occurs. Subsequently, the drop tendency is relatively much milder for FD = 20 mm–50 mm. For FD = 30 mm and θ = 10°–30°, the increase of the maximal normalized acoustic intensities around the skull without mode conversions is faster than that with mode conversions. For θ = 10°–20°, the maximal normalized acoustic intensities around the skull with mode conversions are still higher than those without mode conversions, and both remain lower than the focusing gains. However, when θ = 30°, the maximal normalized acoustic intensity around the skull without mode conversions exceeds the one with mode conversions and both are stronger than the corresponding focusing gains; the reason is there is more reflected acoustic energy without mode conversions than that with mode conversions.

Figure 8(d) shows the ratios of the maximal normalized acoustic intensity around the skull to the focusing gain versus computational model number. Except for θ = 30°, the ratios with mode conversions are greater than those without mode conversions and both the ratios for the other computational models are all smaller than 1. For θ = 0° and FD = 10 mm–50 mm, the ratios decrease, while for FD = 30 mm and θ = 10°–30°, the ratios increase. Specially, for θ = 30°, the ratio without mode conversions is 2.9, which is higher than 2.66 with mode conversions. When θ increases to a larger value, it is believed that the ratio without mode conversions will further exceed the ratio with mode conversions, since the focusing gain with mode conversions will exceed that without mode conversions more, while the maximal normalized acoustic intensity around the skull without mode conversions will exceed that with mode conversions more as analyzed above.

For the 256-element phased array with a uniform emitting acoustic pressure amplitude, numerical results of the normalized acoustic intensities above indicate that a focus designed to be closer to the skull or a larger central incident angle from the array to the skull make it easier to cause serious acoustic energy to deposit around the skull, which is detrimental to the tcFUS therapy. For a central incident angle smaller than 30°, an overestimation of acoustic energy deposition around the focus and an underestimation of acoustic energy deposition around the surfaces of the skull will simultaneously occur when mode conversions are neglected. Nevertheless, when the central incident angle exceeds 30°, the situation may be reversed. We suggest that a location deep in the brain combined with a central incident angle close to 0° may be more suitable for ultrasound emitted from the array to be focused at.

The duration of the tcFUS field is set to be TD = 1 s and the emitting pressure is chosen to be . The temperature rise distributions without and with mode conversions on the cross sections in Fig. 3 are shown in Figs. 9 and 10, respectively. The temperature rise distributions without and with mode conversions on the central axis are acquired from Figs. 9 and 10 for comparison as shown in Fig. 11.

Fig. 9.

	Figure Option ViewDownloadNew Window
	Fig. 9. (color online) Temperature rise distributions ( ) on the cross section without mode conversions, where ((a)–(e)) θ = 0° and FD increases from 10 mm to 50 mm in steps of 10 mm; and ((f)–(h)) FD = 30 mm and θ increases from 10° to 30° in steps of 10°.

Fig. 10.

	Figure Option ViewDownloadNew Window
	Fig. 10. (color online) Temperature rise distributions ( ) on the cross section with mode conversions, where ((a)–(e)) θ = 0° and FD increases from 10 mm to 50 mm in steps of 10 mm, and ((f)–(h)) FD = 30 mm and θ increases from 10° to 30° in steps of 10°.

Fig. 11.

	Figure Option ViewDownloadNew Window
	Fig. 11. (color online) Temperature rises on the central axis with (red lines) and without (green lines) mode conversions, where ((a)–(e)) θ = 0° and FD increases from 10 mm to 50 mm in steps of 10 mm, and ((f)–(h)) FD = 30 mm and θ increases from 10° to 30° in steps of 10°.

Like SubSection 3.1, the temperature rise at the focus and the maximal temperature rise around the skull for each of all computational models are picked up from Figs. 9–11, and plotted in Figs. 12(a) and 12(b), respectively. The drop percentages of the temperature rises at the focuses with mode conversions compared with those without mode conversions are given in Fig. 12(c), whilst the ratios of the maximal temperature rises around the skull to the temperature rises at the focuses are shown in Fig. 12(d). Regardless of the specific values, the variation rules of the curves in Fig. 12 are substantially coincident with the variation rules of the corresponding rules in Fig. 8, where NAI corresponds to TR. The ultrasound field maps well to the ultrasound-induced temperature field. Therefore, in the following, the descriptions related with the variation rules of the curves in Fig. 12 for the temperature rise are not exhaustively repeated, whilst the distinctions between Figs. 8 and 12 are depicted.

Fig. 12.

Figure Option ViewDownloadNew Window

Fig. 12. (color online) Y: with mode conversions; N: without mode conversions; TR: temperature rise. Computational model numbers 1–8 represent the computational models in Figs. 3(a)–3(h), respectively. (a) Temperature rise at the focus, (b) maximal temperature rise around skull, (c) drop percentage of the temperature rise at focus with mode conversions compared with one without mode conversions, and (d) ratios of the maximal temperature rise around the skull to temperature rise at focus.

First of all, figures 9–11 clearly give an intuitive impression that the temperature rise in the skull is higher than that in the focal region in either computational model, since the absorption coefficient of the skull is much greater than those of the water and the brain as given in Table 1.

Figure 12(c) indicates that the average drop percentage of the temperature rises at the focuses with mode conversions compared with those without mode conversions for all computational models is 12.02%. For θ = 0° and FD = 10 mm–50 mm, the temperature rises at the focuses with and without mode conversions, which decline slowly, are respectively around 0.58 °C and 0.65 °C. The average drop percentage is 10.67% as shown in Fig. 12(c). For FD = 30 mm and θ = 10°–30°, the temperature rise at the focus comes down rapidly from 0.6 °C to 0.16 °C without mode conversions and from 0.53 °C to 0.13 °C with mode conversions. The average drop percentage is 14.25%.

In Fig. 12(b), of the maximal temperature rises around the skull for all computational models, the smallest ones with and without mode conversions occur when θ = 0° and FD = 50 mm and are respectively 1.29 °C and 1.18 °C, which are higher than any temperature rise at the focus in Fig. 12(a). From Fig. 12(d), we can more clearly observe that the minimal ratios of the maximal temperature rises around the skull to the temperature rises at the focuses with and without mode conversions are 2.46 and 1.94 respectively, when θ = 0° and FD = 50 mm. Both minimal ratios here are much greater than the corresponding minimal ratios of 0.078 and 0.038 in Fig. 8(d). Thus, the large disparity between the absorption coefficient of the skull and those of the water and the brain will aggravate the difference between the temperature rises around the skull and in the focal region and then may deteriorate the skull overheating problem, even though the acoustic intensity around the skull is much lower than the acoustic intensity in the focal region.

Likewise, an overestimation of the temperature rise around the focus and an underestimation of the temperature rise around the surfaces of the skull simultaneously exist for a central incident angle smaller than 30° when mode conversions are neglected. A reversed situation may also occur if the central incident angle is beyond 30°. Nevertheless, comparing with numerical results of the acoustic intensities in SubSection 3.1, the big characteristic or difference for numerical results of the ultrasound-induced temperature rises is that the heat deposition problem is greatly exacerbated due to the strong absorption of the skull to the acoustic energy. Therefore, such a planar phased array with a uniform emitting acoustic pressure amplitude is more applicable to the non-thermal tcFUS therapy in order to avoid the skull heating problem, especially when the focus is deep in the brain and the central axis of the array is perpendicular to the skull as much as possible. As for the thermal tcFUS therapy utilizing this array, it is better that the microbubbles should be used to deliberately enhance the energy deposition at the focus, in order to relieve the skull heating problem.

Equation (14), which associates with the Kelvin–Voigt equation and the Pennes equation, indicates that the acoustic intensity and the temperature rise in the skull contributed by compressional or shear waves can be calculated separately. To further investigate the contributions of shear waves to the acoustic intensity and the temperature rise in the skull for all computational models in Fig. 3, the spatial distributions of the proportions of the acoustic intensities and the temperature rises from shear waves to the corresponding total ones in Figs. 6 and 10 are given in Figs. 13 and 14, respectively. Figure 15(a) shows the maximal proportions of the acoustic intensities and the temperature rises from shear waves in the skull.

Fig. 13.

	Figure Option ViewDownloadNew Window
	Fig. 13. (color online) Spatial distributions of the proportions of the acoustic intensities from shear waves to the total one, where ((a)–(e)) θ = 0°, FD increases from 10 mm to 50 mm in steps of 10 mm, and ((f)–(h)) FD = 30 mm and θ increases from 10° to 30° in steps of 10°.

Fig. 14.

	Figure Option ViewDownloadNew Window
	Fig. 14. (color online) Spatial distributions of the proportions of the temperature rises induced by shear waves to the total one, where ((a)–(e)) θ = 0° and FD increases from 10 mm to 50 mm in steps of 10 mm, and ((f)–(h)) FD = 30 mm and θ increases from 10° to 30° in steps of 10°.

Fig. 15.

	Figure Option ViewDownloadNew Window
	Fig. 15. (color online) Computational model number dependent (a) maximal and (b) average proportions of the acoustic intensity and the temperature rise in the skull caused by shear waves to corresponding total ones.

From Figs. 13 and 14, the regions in the skull where shear waves contribute more acoustic intensity and temperature rise are mainly located near both sides of the skull, as the incident angles of the ultrasound from the array to these regions are relatively large ( ). A larger incident angle at the surface of the skull is easier to induce more ultrasound waves in the water to be converted into shear waves in the skull. Nevertheless, Fig. 14 has some difference from Fig. 13. Because of the heat conduction effect in biological media, the maximal proportions of the temperature rises from shear waves in the skull are lower than the maximal proportions of the acoustic intensities from shear waves in the skull. As seen clearly in Fig. 15(a), the average of the maximal proportions of the temperature rises from shear waves is 61.2%, which is 35.2% smaller than the average of 94.4% for the acoustic intensities. In addition, the bright area of Fig. 14 is larger than that of Fig. 13, meaning that part of the heat in the skull from shear waves is diffused into the water and the brain adjacent to the skull. It is worth noting that the maximal proportion of the temperature rise from shear waves when FD = 30 mm and θ = 30° is 90.5%, which is closer to the maximal proportion of 99% for the acoustic intensity, in comparison with the situations for the other computational models. When θ = 30°, the incident angle of the ultrasound from the array to the upper side of the skull is greater than 30°. Thus, the region with large proportions ( ) of the acoustic intensities from shear waves in Fig. 13(h) is much greater than those in the other computational models, leading to a relatively slow heat conduction process in each of these regions. Therefore, the maximal proportions of the acoustic intensity and the temperature rise from shear waves for θ = 30° are closer to each other.

Moreover, figure 13 gives us an overall impression that the acoustic intensity in the skull from shear waves is lower than that from compressional waves, since the proportion of the acoustic intensity from shear waves in most of the region of the skull in either computational model is smaller than 50%, and so does Fig. 14 about the temperature rise. Quantitatively, based on Figs. 13 and 14, the average proportions of the acoustic intensities and the temperature rises from shear waves in the skull are calculated from

According to Fig. 15(b), the average proportions of the acoustic intensities and the temperature rises in the skull from shear waves are almost the same for either computational model in Fig. 3, only with a tiny difference. For all the computational models, shear waves averagely contribute 20.54% of the total acoustic intensity and 20.74% of the total temperature in the skull. For θ = 0° and FD = 10 mm–50 mm, the average proportion of the acoustic intensity in the skull from shear waves fluctuates between 16% and 22.06% with an average value of 19.2%, whilst the one of the temperature rise changes between 16.9% and 22.05% with an average value of 19.4%. For FD = 30 mm and θ = 10°–30°, the average proportion of the acoustic intensity in the skull from shear waves increases from 20.65% to 24.94% with an average value of 22.7%, and that of the temperature rise increases from 20.53% to 24.98% with an average value of 22.9%. Both average values for FD = 30 mm and θ = 10°–30° are higher than the corresponding ones for θ = 0° and FD = 10 mm–50 mm, which further indicates that more ultrasound in the water will be converted into shear waves in the skull when the central incident angle increases.

Furthermore, as just mentioned above, the difference between the average proportions of the acoustic intensity and the temperature rise in the skull from shear waves is indeed existent yet tiny and inconspicuous, as shown in Fig. 15(b). However, the empirical formulas under Eq. (17) indicate that the absorption of the skull to shear waves is stronger than that to compressional waves. It seems that the average proportion of the temperature rise in the skull from shear waves should be a little higher than that of the acoustic intensity. We consider that the duration of the ultrasonic field (namely TD) is the key. Initially, except in the skull, there is no acoustic energy from shear waves in the water nor in the brain, while the acoustic energy in the water and the brain is entirely from compressional waves. Thus, the gradient between the energy in the skull and that in the water and the brain from shear waves is larger than that from compressional waves, leading to the heat in the skull from shear waves to diffuse into the water and the brain more rapidly than the heat in the skull from compressional waves. We deem that the average proportion of the temperature rise in the skull from shear waves could be greater than, or equal to, or less than that of the acoustic intensity, with the increase of TD. Perhaps TD = 1 s happens to be used above, and thus leading the average proportions of the acoustic intensity and the temperature rise to be close to each other.

Taking the computational model of FD = 30 mm and θ = 0° for example, we give the spatial distributions of the proportions of the temperature rises induced by shear waves with TD = 0.01, 0.5, 2.0, and 3.0 s in Fig. 16. Figure 16(a) is similar to the corresponding spatial distributions of the proportions of the acoustic intensities from shear waves in Fig. 13(c), since TD = 0.01 s is so short that the heat conduction process has only just begun and is insufficient. With the extension of TD, the heat in the skull from shear waves gradually diffuses into the water and the brain and the maximal proportion of the temperature rise in the skull caused by shear waves drops rapidly from 98.1% to 56.5%, as can be clearly observed in Figs. 16 and 17(a), respectively. For the average proportions of the temperature rise in the skull caused by shear waves, Figure 17(b) shows that it decreases slowly from 22.8% to 21.27% as TD increases from 0.01 s to 3.0 s, and that the average proportion 22.06% for the corresponding acoustic intensity is between 22.8% and 21.27%. Therefore, it is demonstrated that the average proportion of the temperature rise in the skull from shear waves is indeed greater than, equal to and less than that of the acoustic intensity as TD gradually increases.

Fig. 16.

	Figure Option ViewDownloadNew Window
	Fig. 16. (color online) FD = 30 mm and θ = 0°. Spatial distributions of the proportions of the temperature rises ( ) induced by shear waves to the total one, where (a) TD = 0.01 s, (b) TD = 0.5 s, (c) TD = 2.0 s, and (d) TD = 3.0 s.

Fig. 17.

	Figure Option ViewDownloadNew Window
	Fig. 17. (color online) TD-dependent (a) maximal and (b) average proportions of the temperature rise in the skull caused by shear waves to corresponding total ones.